home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / hlvector.lha / hl_vector / ali.cc next >
C/C++ Source or Header  |  1993-08-08  |  8KB  |  277 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *              Numerical Math Package
  6.  *             Aitken-Lagrange interpolation
  7.  *
  8.  * The package allows one to interpolate function value for a given argument
  9.  * value using function values tabulated over either uniform or non-uniform
  10.  * grid. The latter is specified by the vector of node point abscissae.
  11.  * The uniform grid is specified by the grid mesh and the abscissa of
  12.  * the first grid point.
  13.  *
  14.  * Synopsis
  15.  *    double ali(q,x0,s,y)
  16.  *    double q             Argument value specified
  17.  *    double x0             Abscissae for the 1. grid point
  18.  *    double s             Grid mesh, >0
  19.  *    VECTOR y            Vector of function values
  20.  *                    tabulated at points
  21.  *                    x0 + s*(i-y.q_lwb()))
  22.  *                                      The vector has to contain at 
  23.  *                    least 2 elements
  24.  *
  25.  *    double ali(q,x,y)
  26.  *    const double q            Argument value specified
  27.  *    const VECTOR x            Vector of grid node abscissae
  28.  *    const VECTOR y            Vector of function values
  29.  *                    tabulated at points x[i]
  30.  *                                      The vector has to contain at 
  31.  *                    least 2 elements
  32.  * Output
  33.  *    Both functions return the interpolated function value at point q
  34.  *    Interpolation is terminated either
  35.  *        - if the difference between two successive interpolated
  36.  *          values is absolutely less than EPSILON
  37.  *        - if the absolute value of this difference stops
  38.  *          diminishing
  39.  *        - after (N-1) steps, N being the no. of elements in vector y
  40.  *
  41.  * Algorithm
  42.  *    Aitken scheme of Lagrange interpolation
  43.  *    Algorithm is described in the book
  44.  *        "Mathematical software for computers", Institute of
  45.  *         Mathematics of the Belorussian Academy of Sciences,
  46.  *         Minsk, 1974, issue #4,
  47.  *         p. 146 (description of ALI, DALI subroutines)
  48.  *         p. 180 (description of ATSE, DATSE subroutines)
  49.  *    The book essentially describes the IBM SSP package.
  50.  *    
  51.  ************************************************************************
  52.  */
  53.  
  54.  
  55. #include "LinAlg.h"
  56. #include "math_num.h"
  57. #include <std.h>
  58.  
  59. #pragma implementation
  60.  
  61.  
  62. /*
  63.  *------------------------------------------------------------------------
  64.  *             Class that handles the interpolation
  65.  */
  66.  
  67. static class ALInterp
  68. {
  69.   Vector arg(1);        // [1:n] Arranged table of arguments
  70.   Vector val(1);        // [1:n] Arranged table of function values
  71.   double q;                   // Argument value the function is to be
  72.                 // interpolated at
  73.  
  74. public:
  75.                 // Construct the arranged tables for the
  76.                 // uniform grid
  77.   ALInterp(const double q, const double x0, const double s, const Vector& y);
  78.                 // Construct the arranged tables for the
  79.                 // non-uniform grid
  80.   ALInterp(const double q, const Vector& x, const Vector& y);
  81.  
  82.   ~ALInterp(void)    {}
  83.  
  84.   double interpolate(void);    // Perform actual interpolation
  85. };
  86.  
  87. /*
  88.  *------------------------------------------------------------------------
  89.  *
  90.  *        Arranging data for the Aitken-Lagrange interpolation
  91.  *
  92.  * Abscissae (arg) and ordinates (val) of the grid points should be enumerated
  93.  * so that the distance abs(q-arg[i]) would increase as i increases. 
  94.  * Here q is the point the function is to be interpolated at.
  95.  *
  96.  */
  97.  
  98.                 // Construct the arranged tables for the
  99.                 // uniform grid
  100. ALInterp::ALInterp(const double q, const double x0,
  101.            const double s, const Vector& y)
  102. {
  103.   const int n = y.q_no_elems();
  104.   assure( n > 1, "Vector y (function values) has to have at least 2 points");
  105.   assure( s > 0, "The grid mesh has to be positive");
  106.  
  107.   arg.resize_to(n);
  108.   val.resize_to(n);
  109.   ALInterp::q = q;
  110.  
  111.             // First find the index of the grid node which
  112.             // is closest to q. Assign index 1 to this
  113.             // node. Then look at neighboring grid nodes
  114.             // and assign indices to them
  115.             // (kind of breadth-first search)
  116.   int js = (int)( (q-x0)/s + 1.5 );    // Index j for the point x0+s*j
  117.                     // which is the closest to q
  118.   if( js < 1 )                // Check for the case of extrapolation
  119.     js = 1;                // to the left end
  120.   else if( js > n )
  121.     js = n;                // or to the right end
  122.  
  123.   register int dir;            // Direction to the next closest
  124.                     // to q grid node
  125.   dir = ( q > x0 + (js-1)*s ? 1 : -1 );
  126.  
  127.   register int jcurr = js, jleft = js, jright = js;
  128.   register int i;
  129.   for(i=1; i<=n; ++i)            // Pick up elements x0+s*i
  130.   {                    // in the neighborhood of q
  131.     arg(i) = x0 + (jcurr-1)*s;
  132.     val(i) = y(jcurr-1+y.q_lwb());    // Once the closest to q point js
  133.     if( jright >= n )            // is found, we pick up points
  134.       dir = -1;                // alternatively to the right
  135.     if( jleft <= 1 )            // and to the left of the js
  136.       dir = 1;                // further and further
  137.     if( dir > 0 )
  138.     {
  139.       jcurr = ++jright;
  140.       dir = -1;
  141.     }
  142.     else
  143.     {
  144.       jcurr = --jleft;
  145.       dir = 1;
  146.     }
  147.   }
  148. }
  149.  
  150.  
  151. /*
  152.  * The function that defines the order while sorting the array of indices
  153.  * in array arg using qsort
  154.  * The function is given two pointers to elements in indices array, i.e.
  155.  * to the two indices for the _XtoQ array, _XtoQ[i] being abs(x[i]-q).
  156.  * The function returns -1, 0, or +1 depending on the fact if 
  157.  * abs(x[i]-q) is less, equal, or greater than abs(x[j]-q)
  158.  */
  159.  
  160. static Vector * _XtoQ;
  161.  
  162. int compare_indices(const short *ip, const short *jp)
  163. {
  164.   register double delta = (*_XtoQ)(*ip) - (*_XtoQ)(*jp);
  165.   if( delta < 0 )
  166.     return -1;
  167.   else if(delta == 0)
  168.     return 0;
  169.   else
  170.     return 1;
  171. }
  172.  
  173.                 // Construct the arranged tables for the
  174.                 // non-uniform grid
  175. ALInterp::ALInterp(const double q, const Vector& x, const Vector& y)
  176. {
  177.   const int n = y.q_no_elems();
  178.   assure( n > 1, "Vector y (function values) has to have at least 2 points");
  179.   are_compatible(x,y);
  180.  
  181.   arg.resize_to(n);
  182.   val.resize_to(n);
  183.   ALInterp::q = q;
  184.                     // Selection is done by sorting x,y arrays
  185.             // in the way mentioned above. Fisrt an array
  186.                     // of indices is created and sorted, then arg,
  187.             // val arrays are filled in using the sorted indices
  188.    
  189.   short indices[n];
  190.   Vector xtoq(0,n-1);
  191.   register int i;
  192.  
  193.   for(i=0; i<n; i++)            // 0..n-1 correspond to the
  194.   {                    // unsorted x
  195.     indices[i] = i;
  196.     xtoq(i) = fabs(q - x(i+x.q_lwb()));
  197.   }
  198.   _XtoQ = &xtoq;
  199.                         // Sort indices so that
  200.                     // |q-x[ind[i]]| < |q-x[ind[j]]|
  201.                     // for all j>i
  202.   qsort(indices,n,sizeof(indices[0]),compare_indices);
  203.  
  204.   for(i=arg.q_lwb(); i<=arg.q_upb(); i++)
  205.   {
  206.     register int ind = indices[i-arg.q_lwb()];
  207.     arg(i) = x(x.q_lwb() + ind);
  208.     val(i) = y(y.q_lwb() + ind);
  209.   }
  210. }
  211.  
  212. /*
  213.  *------------------------------------------------------------------------
  214.  *            Aitken - Lagrange process
  215.  *
  216.  *  arg and val tables are assumed to be arranged in the proper way
  217.  *
  218.  */
  219.  
  220. double ALInterp::interpolate()
  221. {
  222.   register double valp = val(1);    // The best approximation found so far
  223.   register double diffp = MAXDOUBLE;    // abs(valp - prev. to valp)
  224.   register int i,j;
  225.  
  226. #ifdef DEBUG
  227.   arg.print("arg - interpolation nodes");
  228.   val.print("Arranged table of function values");
  229. #endif
  230.             // Compute the j-th row of the Aitken scheme and
  231.             // place it in the 'val' array
  232.   for(j=2; j<=val.q_upb(); j++)
  233.   {
  234.     register double argj = arg(j);
  235.     register REAL&  valj = val(j);
  236.     for(i=1; i<=j-1; i++)
  237.       valj = ( val(i)*(q-argj) - valj*(q-arg(i)) ) / (arg(i) - argj);
  238.  
  239. #ifdef DEBUG
  240.     message("\nval(j) = %g, valp = %g, arg(j) = %g",valj,valp,argj);
  241. #endif
  242.  
  243.     register double diff = fabs( valj - valp );
  244.  
  245.     if( j>2 && diff == 0 )          // An exact result has been achieved
  246.       break;
  247.  
  248.     if( j>4 && diff > diffp )        // Difference stops diminishing
  249.       break;                // after the 4. step
  250.  
  251.     valp = valj;  diffp = diff;
  252.   }
  253.  
  254.   return valp;
  255. }
  256.  
  257.  
  258.  
  259. /* 
  260.  *=======================================================================
  261.  *                Root modules
  262.  */
  263.  
  264.                 // Uniform mesh x[i] = x0 + s*(i-y.lwb)
  265. double ali(const double q, const double x0, const double s, const Vector& y)
  266. {
  267.   ALInterp al(q,x0,s,y);
  268.   return al.interpolate();
  269. }
  270.  
  271.                 // Nonuniform grid with nodes in x[i]
  272. double ali(const double q, const Vector& x, const Vector& y)
  273. {
  274.   ALInterp al(q,x,y);
  275.   return al.interpolate();
  276. }
  277.